Universality-class dependence of energy distributions in spin glasses 
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We study the probability distribution function of the ground-state energies of the disordered one- 
dimensional Ising spin chain with power-law interactions using a combination of parallel tempering 
Monte Carlo and branch, cut, and price algorithms. By tuning the exponent of the power-law 
interactions we are able to scan several universality classes. Our results suggest that mean-field 
models have a non-Gaussian limiting distribution of the ground-state energies, whereas non-mean- 
field models have a Gaussian limiting distribution. We compare the results of the disordered one- 
dimensional Ising chain to results for a disordered two-leg ladder, for which large system sizes can 
be studied, and find a qualitative agreement between the disordered one-dimensional Ising chain in 
the short-range universality class and the disordered two-leg ladder. We show that the mean and 
the standard deviation of the ground-state energy distributions scale with a power of the system 
size. In the mean-field universality class the skewness does not follow a power-law behavior and 
converges to a nonzero constant value. The data for the Sherrington-Kirkpatrick model seem to 
be acceptably well fitted by a modified Gumbel distribution. Finally, we discuss the distribution 
of the internal energy of the Sherrington-Kirkpatrick model at finite temperatures and show that 
it behaves similar to the ground-state energy of the system if the temperature is smaller than the 
critical temperature. 

PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q 



I. INTRODUCTION 



Averages of physical quantities and their fluctuations 
play an important role in statistical physics; however the 
knowledge of the "average" behavior of a quantity often 
does not provide sufficient information to fully charac- 
terize a system, especially if the probability distribution 
of the quantity in question is non-Gaussian, e.g., when it 
has a nonvanishing skewness. Hallmark examples of such 
distributions are power-law or exponential distributions, 
which in nature occur in relation to earthquakes^ mag- 
netic fluctuations) 2 ' stock markets, 3 directed polymers in 
a random medium^*^ coauthorships in publications, the 
Internet, and other complex networks.? Many of these 
systems are characterized by the absence of a character- 
istic length scale such that rare events involving large 
parts of the system become important and strongly in- 
fluence the average of various quantities. 

Recently, there has been a renewed interest in the 
ground-state energy distribution P{E) and its limiting 
form P oa {E) of the mean-field Sherrington-Kirkpatrick 
(SK) spin-glass modeffi^iSi and of short-range spin 
glasses in two and three dimensions^ While studies of 
the mean-field model have found a non-Gaussian limiting 
distribution^^ the study of small system sizes of two- 
and three-dimensional short-range spin glasses^ have 
found a Gaussian limiting distribution in the thermo- 
dynamic limit. This is supported by the fact that sys- 
tems with short-range interactions can be subdivided into 
smaller subsystems, coupled weakly enough to contribute 
almost independently to the total energy and leading to 



a Gaussian distribution via the central-limit theorem; 11 
however it is important to note that the weak coupling 
between the ground-state energies of subsystems below 
or at an ordering temperature is not self-evident. 

Our goal is to consolidate the different limiting cases of 
short-range and long-range interactions in spin glasses 12 
by studying a disordered one-dimensional Ising spin chain 
with power-law interactionsii 3 ^^^^^ The model has 
the advantage over conventional models in that by tun- 
ing the power-law exponent, several universality classes 
ranging from mean-field type behavior to a short-range 
spin glass can be probed for a large range of system 
sizes. We show that the presence or absence of mean- 
field behavior 1 - is reflected in the limiting distribution of 
the ground-state energies. We also study a two-leg short- 
range spin ladder, where an exact transfer-matrix algo- 
rithm can be applied, in order to compute the ground- 
state energy distribution for large system sizes and to 
obtain a comparison for the results of the disordered 
Ising chain with power-law interactions in the short-range 
phase. Using a large range of system sizes, our results 
clearly show that mean-field spin-glass models have a 
non-Gaussian limiting distribution with a finite skewness 
in the thermodynamic limit, whereas the limiting dis- 
tributions for nonmean-field models are Gaussian (also 
referred to as "Normal" ) . In addition, we also find that 
the distribution of the internal energy of mean-field mod- 
els is non-Gaussian if the temperature is lower than the 
critical temperature. 

We do not attempt to make a prediction regarding 
the exact functional form of the limiting distribution for 
an arbitrary spin-glass model. Bouchaud et alM> have 
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shown for small system sizes that typical short-range 
spin-glass models have a Normal limiting distribution. 
This is not the case for the mean-field model, thus pos- 
ing the question of whether the limiting distribution falls 
into one of the standard three universality classes for the 
minimum of uncorrelated variables: 19 Gumbel, Fischer- 
Tippet-Frechet, and Weibull distributions. The results of 
Bouchaud et almL cannot determine with certainty which 
limiting distribution fits the data bestiSi Our results sug- 
gest that a modified Gumbel distribution2iL2I fits the 
data for the SK model best, although a detailed probing 
of the tails of the energy distributions would be required 
to make a definite statement if corrections to the modi- 
fied Gumbel distribution are required. For finite values of 
the power-law exponent we add a quadratic correction to 
the modified Gumbel distribution and show that for finite 
system sizes the data are well described by this function. 
In addition, the quadratic correction (Gaussian) domi- 
nates for increasing system size in the nonmean-field uni- 
versality class, thus showing that in the thermodynamic 
limit a Normal distribution is recovered. 

General scaling arguments are presented in Sec. [H] In 
Sec. IIIII we present results on a one-dimensional two- leg 
ladder with short-range interactions in order to illustrate 
the expected results for a short range model for very 
large system sizes. Results on the one-dimensional Ising 
spin chain with power-law interactions at zero and finite 
temperatures are presented in Sec. II VI We conclude in 
Sec. The numerical methods used to compute the 
ground-state energies22iS are described in the Appen- 
dices. 



II. STATISTICAL DESCRIPTION OF DATA 



The skewness of a distribution of M values {Ei} is given 
by 
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where (E) and oe are given by 
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respectively. Note that the skewness is a dimensionless 
quantity. Following previous results by Ref . ^] we expect 
the skewness to decay as 



Ce = ci + c 2 L~ 



(6) 



with 7 > 0. As we shall see later, ci = for the short- 
range limit of the model. We also want to test whether 
the scaled probability distribution functions P(e) with 
e = {E — (E))/(Te converge to a limiting form Poo(e) for 
L — > oo. If this is the case, then data for the ground-state 
energies should be scalable via 



P(E) 



1 ( E-{E) 



(7) 



where (E) and oe are given by Eqs. Q and (JSJ, respec- 
tively. 



In general, we expect the ground-state energy of a dis- 
ordered system to be a random variable with mean (E), 
standard deviation erg, and skewness Ce^* I n this work 
we study the size dependence of the aforementioned ob- 
servables. In particular, we make the ansatz that the 
mean ground-state energy of a (one-dimensional) random 
system scales as 

(E)/L = eoc +aL~" , (1) 

where L represents the system size (and number of spins) 
and uj describes the leading finite-size corrections for the 
energy per spin. We keep the extra factor of L in Eq. QJ, 
as well as in the following definitions in order to be able 
to compare to the exponent estimates of Ref. 0- The 
standard deviation of the ground-state energy of a gen- 
eral disordered system can be expected to be determined 
by an exponent p via 

ue/L — bL~ p . (2) 



III. TWO-LEG SPIN-GLASS LADDER 

To compare results for the one-dimensional Ising spin 
chain with power-law interactions with a simple bench- 
mark model for which large system sizes can be studied, 
we consider a disordered (short-range) Ising model on a 
two- leg ladder (see Fig. QJ. The couplings J,j between 
nearest-neighbor spins are chosen from a Normal distri- 
bution with zero mean and unit standard deviation. A 
system of length L is described by the Hamiltonian 

L 

T~L = / J J(l,a),(l,b) S(l,a) S(l,b) 
1=1 
L-l 

+ X] J (l,i),(l+hi) S (l+hi)' ( 8 ) 

1=1 i=a.b 

The first summation in Eq. JSJ runs over all rungs Z, while 
the second summation runs over all exchanges between 
the rungs, and Sn^ = ±1 is the value of the (Ising) 
spin on the ith leg of the Ith rung of the ladder. The 
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FIG. 1: Illustration of a step in the transfer-matrix calcula- 
tion. Starting with a system of size L (panel I) whose ground- 
state energy is known as a function of the spins (open circles), 
we add another rung and calculate the change in energy Ae 
given by the dashed line [see panel II and Eq. Ulul l. The 
ground-state energy of the system as a function of the spins 
in the L+ l-th rung is then calculated by taking the minimum 
of Eg + Ae over all configurations of the Lth rung [see panel 
III and Eq. ©]. 



ground-state energy of the system can be efficiently calcu- 
lated with a transfer-matrix algorithmi2SiS& The transfer- 
matrix algorithm computes the ground-state energy of a 
system of size L in 0{L) time so that large systems can 
be studied. The disorder average has to be performed 
explicitly by repeating the algorithm for a number of dis- 
order realizations. 



A. Numerical Method: Transfer Matrices 

We can explain the transfer matrix algorithm by 
starting with a ladder of length L and assuming that 
the ground-state energy E g (L, {S(l^}) of the ladder is 
known as a function of the spin configuration {S( L ^} of 
the Lth rung. We add the spins of the L + 1-th rung to 
the system as illustrated in Fig. ^ and use the relation 

E K (L + 1,{S (L+1A }) = min [E K (L,{S (LA }) 

+ A E ({S (Lti) },{S iL+1>i) })](9) 

to integrate out the spins of the L-th rung and to obtain 
Eg as a function of the spins of the L + l-th rung. Here 

^(1%,!)}, {S(L + l,i)}) = /J(L,i),(L+l,i) S(L,i) S(L+l,i) 

i—a,b 

+ 'J '(L+l, a), {L + l, b) S (L+l,b) (10) 

is the exchange energy of the spins added on the L + 1- 
th rung with themselves and with the spins of the Lth 
rung. Starting with two spins, we iterate this procedure 
until the system has the desired size L max . The final 
ground-state energy is then obtained by minimizing over 



the spins of the last rung 

E g (L max ) = min [E s (L maX) {£(z, max ,i)})] ■ (H) 

We repeat the calculation until a desired number of dis- 
order realizations is obtained. 



B. Results 

In Fig. we scale the data for the energy of the lad- 
der system according to Eq. Q for system sizes up to 
L = 10 4 . For each system size we compute 10 6 samples. 
The data scale well, although deviations are present in 
the tails. In particular, for small L the distribution is 
clearly skewed. For the short-ranged ladder system we 
obtain a clear power-law decay of the skewness accord- 
ing to Eq. Jfjjl with 7 w 0.5 (and c\ = 0), as can be seen 
in Fig. |3J This suggests that in the thermodynamic limit 
the ground-state energies are Gaussian distributed. For 
completeness, we quote the results for the size depen- 
dence of the mean and standard deviation. We obtain 
for the mean energy 



(E)/L = -2.12582(8) - 0.801(8)L 



and thus DRil. For the fluctuations 



-0.996(4) 



o E = 0.976(5)L ( 



0.497(8) 



(12) 



(13) 



Our results therefore show that p ~ —1/2 as in the case of 
the one-dimensional Ising chain (see below) , and (E) / L — 
eoo ~ 1/L. 

The scaling of the skewness to zero with a power law 
and the results for very large system sizes already suggest 
that for short-ranged systems the limiting distribution 
of the ground-state energy is Normal. This result dif- 
fers from recent results^ for the mean-field Sherrington- 
Kirkpatrick model^S 7 . where the limiting distribution 
seems to have a finite skewness and thus cannot be prop- 
erly described by a Gaussian. Hence, it is desirable to 
study a system that allows to interpolate between both 
cases to verify whether the change of the distribution 
coincides with a general change of the universality class. 
This is indeed the case for the one-dimensional long-range 
Ising spin glass, which is studied in the following section. 



IV. ID ISING CHAIN 

The Hamiltonian for the one-dimensional long-range 
Ising spin glass with power-law interactions is given by 



H — ^2 JijSiSj , 



(14) 



where Si = ±1 represent Ising spins evenly distributed 
on a ring of length L in order to ensure periodic boundary 
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FIG. 2: (Color online) Scaling of the ground-state energy 
according to Eq. £J f° r the ladder system. The data scale 
well, although deviations in the tails suggest that the skewness 
of the function is changing with system size L. 
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FIG. 3: (Color online) Skewness of the energy distributions 
as a function of system size L for the ladder system. The 
skewness can be well fitted to a power-law decay with an 
exponent 7 ~ 1/2. This suggests that in the thermodynamic 
limit the limiting distribution is Gaussian (zero skewness). 



conditions. The sum is over all spins on the chain and 
the couplings Jy are given bjii& 



(15) 



where the e<y are chosen according to a Gaussian distri- 
bution with zero mean and standard deviation unity 



1 



exp(-e?./2) 



(16) 



and Tij = (L/ir) sm[(w\i — j\) / L] represents the geometric 
distance between the spins on the ring^ The power-law 
exponent a determines the range of the interactions and 
thus the universality class of the model, as described in 
the next section. The constant c(cr) in Eq. I|15|) is chosen 
to give a mean-field (MF) transition temperature T C MF = 
1, where 



T, 



MF 



j^i, i fixed 



ij]av 



i. i fixed l 3 



(17) 



Here [•••]av denotes an average over disorder. In 
this work we compute unsealed energies for the one- 
dimensional Ising chain. Thus we find the optimal con- 
figuration of spins {Si} that minimizes the Hamiltonian 
in Eq. I|14[) for a given set of interactions {Jij}, i.e., 



E = minH({J ij },{S i }). 



(18) 



The (commonly used) energy per degree of freedom e is 
then given by e = E/L. 



A. Phase Diagram 

The d-dimensional long-range Ising spin glass with 
power-law interactions has a very rich phase diagram in 
the d-a plane. This is summarized in Fig. ^ which is 
based on work performed by Bray et a/.— and by Fisher 
and Husei^ who present a detailed analysis of the role of 
long-range interactions within the droplet model. Spin- 
glass behavior is controlled by the long-range part of the 
interaction if a is sufficiently small, and by the short- 
range part if a is sufficiently large. More precisely, one 
has long-range behavior if the stiffness exponent of the 
long-range (LR) universality class #lr is greater than 
that of the short-range (SR) universality class #sr and 
vice versa. In addition, there is an exact result for #lr, 
namely(iii4 



#lr = d-a, 
so long-range behavior occurs if 

a < a c (d) =d- 9 S n{d) 



(19) 



(20) 



Equation l|19|l indicates that critical exponents depend 
continuously on a in the long-range region, even though 
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d i d =2a o c (d) 




11/2 12 a 

FIG. 4: (Color online) Sketch of the phase diagram in the d-a 
plane for the spin-glass state of the disordered long-range Ising 
model with power-law interactions following Ref. |l4. The 
light shaded region (LR + ) is where there is both a finite T c 
and the spin-glass state is controlled by the long-range part 
of the interaction. The thick solid line separates the region of 
short-range behavior (SR) from that of long-range behavior 
and is denoted by a c (d). The thick dashed line separates re- 
gions where T c = (e.g., LR°) from regions where T c > 0, i.e., 
it corresponds to a zero stiffness exponent. The dark shaded 
region (MF, a < d/2) is where there is no thermodynamic 
limit unless the infinite-range interactions are scaled appro- 
priately by the system size. The calculations are performed 
for d = 1 (marked by a horizontal dashed red line), for which 
<7 c (d) = 2 within a droplet picture approximation. These val- 
ues of a are marked. Note that we refer to the infinite-range 
region in the phase diagram as "mean-field region" in order to 
be consistent with previous studies, even though the mean- 
field region extends to d — (2/3)o\ (Figure adapted from 
Ref. El) 



they are independent of a in the region controlled by the 
short-range part of the interaction. Thus we expect to be 
able to tune the different universality classes by changing 
the exponent a. The condition for a finite-temperature 
transition is 9 > 0, where 8 refers here to the greater of 
6?sr and #lr- For the short-range model, there is a finite- 
temperature transition (i.e., #sr > 0) for d larger than 
the lower critical dimension di, which is found numeri- 
cally to lie between 2 and 3i 30 i 31 i 32 i 33 i 34 i 35 i 36 For d = 1, as 
in the present study, we obtain a finite transition temper- 
ature for a < 1. For a < d/2, the model would not have 
a thermodynamic limit (T c would diverge) if the interac- 
tions were not scaled as shown in Eq. 1|17[) . The scaling 
leads to a power-law dependence on L with a negative 
exponent, i.e., c(cr) — > for L — > oo. a — corresponds 
to the SK model and leads to c(0) ~ 1/L. 
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FIG. 5: (Color online) Unsealed energy distributions for sev- 
eral system sizes for the one-dimensional Ising chain with 
a = 0.75 (LR+ phase). 



B. Numerical Methods 

Ground-state energies for the one-dimensional Ising 
chain are computed using the parallel tempering Monte 
Carlo methodiSiSLSS^S when the power-law exponent 
a is small, and the branch, cut, and price (BCP) 
algorithm^fl^i^ when a is large. As reported in Ref. fl7L 
the time to compute a ground-state instance using the 
parallel tempering Monte Carlo method scales in prac- 
tice with a power of the system size for a < 1.25, whereas 
for large values of a the time to compute a ground state 
scales ~ exp(aL), with a a constant. In this case we 
use the BCP algorithm which performs best for short- 
range interactions, thus ideally complementing the par- 
allel tempering method. Details about the algorithms 
used and simulation parameters can be found in the Ap- 
pendices. 



C. Results 

For each system size we compute 10 5 ground-state re- 
alizations for system sizes up to L = 192 (see Table for 
details). In Fig. we show a representative set of the 
unsealed data for a = 0.75 in the LR + -phase for several 
system sizes. Data for other values of a show a similar 
qualitative behavior. The data in Fig. [5] can be scaled 
according to Eq. (JJJ, the result is displayed in Fig. 
The data are clearly skewed and the tails indicate that 
the skewness depends on the system size. In Fig. {7\ we 
show scaled data for the ground-state energy distribu- 
tions for a — (SK limit, MF phase). The data also 
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FIG. 6: (Color online) Scaled ground-state energy distribu- 
tions for several system sizes for the one-dimensional Ising 
chain with a = 0.75 (LR + phase). The dashed vertical line 
is a guide to the eye to illustrate the skewness of the distri- 
bution. The spread of the data in the tails suggests that the 
skewness changes with system size. 



FIG. 7: (Color online) Scaled ground-state energy distribu- 
tions for several system sizes for the one-dimensional Ising 
chain with a = (MF phase, SK model). The dashed ver- 
tical line is a guide to the eye to illustrate the skewness of 
the distribution. The data show little spread in the tails sug- 
gesting a weaker dependence on L than for larger values of 
a. 



show a clear asymmetry, but the spread in the tails is 
noticeably smaller than for larger values of a (see Figs. El 
and |SJ) suggesting a smaller dependence of the skewness 
of the distribution on the system size. 

In order to better quantify the aforementioned behav- 
ior, in Fig. we present data for the skewness as a func- 
tion of system size for several values of the power-law 
exponent a . The data show that for a > 0.5 the skew- 
ness of the ground-state energy distributions decays with 
a power law |Ce| ~ £~ 7 , with 7 w 0.5 in the SR phase, 
whereas for a < 0.5 (MF phase) the skewness is well 
fitted by Eq. JBJl with c\ > thus tending to a con- 
stant in the thermodynamic limit. This means that the 
mean-field models present a singular behavior in which 
the ground-state energy fluctuations are non-Gaussian 
in the thermodynamic limit. This is not the case for the 
nonmean-field universality class where a limiting Gaus- 
sian behavior is obtained for L — ► 00. Note that 7 « 0.5 
for a > 1 for which T c = 0, in agreement with the results 
for the ladder system studied in Sec. IIIII 

We also study the size-dependence of the mean en- 
ergy as a function of a. For the mean-field Sherrington- 
Kirkpatrick models (a = 0) it is known that u> ~ 
2^/3 > 7 i 8 J ii J 43 Q ur resu jts agree well with this prediction, 
i.e., u> = 0.64(1) (the quality of fit probability^ is 
Q = 0.51; the fit is performed for L > 64). Unfortu- 
nately, there are no predictions for the different expo- 
nents for a > 0, thus we will focus on comparing the 
present results to data for the SK model. In Fig. 1101 we 
show data for all values of a studied. For increasing a, to 
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FIG. 8: (Color online) Scaled ground-state energy distribu- 
tions for several system sizes for the one-dimensional Ising 
chain with a = 2.50 (SR phase). The dashed vertical line is a 
guide to the eye to illustrate the skewness of the distribution. 
The data show a moderate dependence on the system size. 
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FIG. 9: (Color online) Skewness £e as a function of system 
size L for several values of a. For a < 0.5 (MF phase) the data 
scale as £g ~ ci + c^L' 1 with c\ > [Eq. JBJ] thus tending 
to a constant in the thermodynamic limit. For a > 0.5 the 
skewness decays with a power-law behavior, i.e., ci = (fits 
done for L > 64). Note that for a > 1, for which T c = 0, 
7 « 0.5 in agreement with the results for the ladder system 
presented in Sec. Mil 



FIG. 10: (Color online) Mean [(E) - E^/L as a function of 
system size L for several values of a. The data are expected 
to decay as a power of the system size with an exponent u. 
Note that for the SK limit w«2/3atT = 0in agreement 
with other predictions (Refs. Fill l7t ITU andOJ [see Fig.lTTlfor 
(jj(ct)] and that for a — 2.5 we obtain u — 2.29(17). 



increases rapidly and then saturates at u> w 2 in the SR 
phase. This can be understood by studying the model for 
a — 00. In this limit there is no frustration, except that 
with a 50 % probability there will be a broken bond due 
to the periodic boundary conditions. Since the weakest 
bond will be broken, for a continuous distribution the en- 
ergy scales as ~ 1/L. Since the total energy scales with 
system size, we expect the finite-size correction to the 
average energy per spin to be ~ l/£ 2 , i.e., ui = 2. This 
behavior can be seen in Fig. llll where we show the behav- 
ior of u> [see Eq. {Q] in detail for the different universality 
classes. 

The behavior of the energy fluctuations is shown in 
Fig. E| as a function of system size L for several values 
of a. In the SK limit there are contradicting predic- 
tions regarding the power-law exponent p of the energy 
fluctuations erg. While Crisanti et al. 44 find p = 5/6, 
Bouchaud et alm^ and Aspelmeier et a/4* find p = 3/4. 
In this work we obtain p = 0.775(2) (Q = 0.58; fits done 
for L > 64), which is also in agreement with the work 
by Palassinii In Fig. ^| we show the a dependence of 
p. It is noteworthy that p decreases from the mean-field 
value ~ 3/4 to 1/2 in the short-range universality class. 
This is to be expected as for a — > 00 the central limit 
theorem predicts that p = 1/2 iAi Note that the results 
found agree with the prediction of the short-range ladder 
system in Sec. IIII Bl 
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FIG. 11: (Color online) Exponent of the mean energy (u>) as 
a function of a, according to Eq. (0. u increases from the 
SK value (~ 2/3) for increasing a. The exponents are only 
estimated for the four largest system sizes studied for a given 
value of <t. See Table for details. In this and following fig- 
ures, the boundaries between the different universality classes 
are denoted by vertical dashed lines. 
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FIG. 12: (Color online) Standard deviation oe as a function 
of system size L for several values of a. The data for oe/L 
are expected to decay as a power of the system size with an 
exponent p. Note that for the SK limit p « 3/4, in agreement 
with other predictions (Ref. EM 0, and IllT) [See Fig. 1131 for 



FIG. 14: (Color online) Difference in area between the ac- 
tual data for the energy probability distributions of the one- 
dimensional Ising chain to a Gaussian limiting distribution as 
a function of system size L for several values of a [see Eq. 121H 1. 
In the MF phase [a < 0.5) the area difference tends to a con- 
stant for increasing system size, whereas for a > 0.5 the area 
difference decays with a power of the system size. Note the 
close resemblance to the behavior found for the skewness of 
the distribution, Fig. [5] 
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FIG. 13: (Color online) Exponent for the energy fluctuations 
p as a function of a [see Eq. ©]■ For a —* p ~ 3/4 in 
agreement with Refs. 14EL llll and0. For a — > oo, p — > 1/2, 
as predicted by the central limit theorem, and in agreement 
with the results on the ladder system presented in Sec. IIII Bl 



D. Limiting Distribution 

In order to further strengthen the conjecture that 
ground-state energy distributions remain skewed in the 
thermodynamic limit for the mean-field phase, in this 
section we study the area deviation of the normalized 
energy distributions in comparison to a Normal distribu- 
tion N(e). We define the area difference A via 

A = JjP(e) - N(e)\de , (21) 

where P(e) are the actual rescaled data [Eq. J7J]. In 
Fig-El we show the area difference as a function of system 
size L for several values of a. The data for a < 0.5 
can be well fitted by a functional form ~ / + g/L h , i.e., 
the area difference tends to a nonzero constant in the 
thermodynamic limit. This is not the case for a > 0.5 
where the area difference decays with a power law of the 
system size, thus showing that the difference between the 
data and a Gaussian limiting distribution decreases for 
increasing L. 

Palassin: 7 has fitted the data for the scaled probability 
distribution functions of the SK model, Eq. (7J), to a mod- 
ified Gumbel distribution2ii2i2i g m (e), and finds good 
agreement between the data and the fit, especially when 
studying the cumulative distributions Q(e) — J e P(x)dx. 
In addition, Palassini shows that the best fit seems to be 
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obtained for m — 6, although to date it is unclear why 
the aforementioned value of m fits the data best. Because 
outside the MF universality class the limiting distribu- 
tion function seems to converge to a Normal distribution, 
we modify the standard modified Gumbel distribution by 
taking into account a Normal contribution^ i.e., 



9' m ( e ) = N (y)9m(y) , 



where 



and 



(22) 



(23) 



flm(e) = wi exp [my - me v ] (24) 
is the modified Gumbel distribution and 
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N(e) = w 2 exp [m 2 y 2 ] (25) 

is a Normal distribution. Here /i is the most-probable 
value, v a standard deviation, and Wi represents an over- 
all normalization factor. For simplicity, we can fix m = 6 
and study the behavior of the coefficient iri2 as a func- 
tion of system size for different values of a. Note that 
for m,2 — g' m {x) oc g rn (x), up to a global scaling fac- 
tor. A multiplicative ansatz [instead of, for example, an 
additive ansatz of the form aN(e) + bg m (e)] can be mo- 
tivated by keeping in mind that for short-range interac- 
tions, the system can be divided into subsystems which 
contribute almost independently to the total energy. In 
Fig. 1151 we show data for ni2 versus L for a few repre- 
sentative values of a. Our results show that for a = 
(SK model) ui2 converges to a value close to zero for 
L — > oo. For < a < 0.5 the limiting distribution is 
non-Gaussian, yet m,2 is small, but finite. For a > 0.5 
the Gaussian contribution via m,2 dominates in the ther- 
modynamic limit (at least for a finite fitting region), as 
can be seen in Fig. ED This shows that the energy dis- 
tributions in the SK model can be well described in the 
thermodynamic limit by a modified Gumbel distribution. 
In order to test the existence of small Gaussian correc- 
tions to the modified Gumbel distribution for the SK 
model, large-scale simulations probing the tails of the 
distribution function in detail would be required which 
are beyond the scope of this work. For all other values 
of a < 0.5 there are clearly Gaussian corrections to the 
Gumbel distribution, whereas for a > 0.5 the data in 
the thermodynamic limit are well described by a Normal 
limiting distribution. This could be due to the fact that 
there are no length scales associated with the mean-field 
model. Thus any length-scale associated effects will scale 
with system size. This is not the case in the short-range 
models where a length scale will not necessarily scale with 
system size, therefore yielding a Normal distribution in 
the thermodynamic limit. 



FIG. 15: (Color online) Coefficient to the quadratic cor- 
rection term in the modified Gumbel distribution, Eq. 1221 . 
as a function of system size for several values of the power- 
law exponent a. The data show that m-2, converges to a value 
close to zero for the SK model in the thermodynamic limit 
thus suggesting that the energy distributions of the SK model 
are possibly well described by a modified Gumbel distribution 
function in agreement with results from Ref. 0. For all a > 0, 
m,2 tends to a finite negative value in the thermodynamic 
limit. For large values of a, U12 dominates thus showing that 
in the SR universality class the limiting probability distribu- 
tion is well described by a Gaussian. The dashed lines are 
guides to the eye. 

E. Finite Temperatures 

We want to test if the fact that the ground-state en- 
ergy distribution of the SK model is skewed in the ther- 
modynamic limit is a unique property of the ground 
state, or if similar effects can be observed at finite tem- 
peratures. Because the parallel tempering Monte Carlo 
method used to compute the ground-state energies of the 
one-dimensional Ising chain at small values of a requires 
the system to be simulated at several temperatures rang- 
ing to values well above the spin-glass transition (for the 
SK model T c — 1), we have also studied the behavior of 
the internal energy distributions in the mean-field limit 
as a function of temperature. The internal energy U for 
a given disorder realization {Jij} is given by 

U = (H({J ij },{S i })) , (26) 

where the Hamiltonian Ti is given by Eq. 1)140. Here (• • • ) 
represents a thermal average over t eq Monte Carlo steps 
that we perform after having equilibrated the system for 
a time t eq (see Table [I] for details). 

Figure El shows data for the skewness of the inter- 
nal energy distributions as a function of system size for 
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FIG. 16: (Color online) Skewness of the internal energy prob- 
ability distribution functions of the SK model as a function 
of system size for different temperatures. The data show a 
curvature for T < T c in a log-log scale thus suggesting that 
the skewness converges to a constant value in the thermody- 
namic limit. For T > T c — 1 the skewness decays with a 
power of the system size (Normal limiting probability distri- 
bution function) . The inset shows the skewness of the internal 
energy distribution of the SK model for L — 192 (largest sys- 
tem size studied) as a function of temperature. The data show 
that for finite system sizes the skewness seems to peak at the 
transition (T c = 1, shaded area). 



several temperatures ranging from the ground-state to 
well above the critical temperature. The results show 
that the skewness of the distributions tend to a constant 
value in the thermodynamic limit for T <T C (curved fit- 
ting functions in a log- log plot, Fig. I16f) . thus showing 
that skewed energy distributions seem to persist for any 
temperature below the critical temperature. For tem- 
peratures above the critical temperature, the skewness 
shows again a power-law behavior thus suggesting that 
for T > T c the limiting distribution is Normal, as one 
would expect. Therefore, the limiting probability dis- 
tribution is skewed in the thermodynamic limit for all 
temperatures below the critical point. 

The inset of Fig. ^| shows the skewness of the proba- 
bility distribution function of the internal energy of the 
SK model for L — 192 as a function of temperature. 
The data show a peak around T c = 1. We expect the 
functional form of the ground-state energy distribution 
to remain approximatively the same for L — > oo when 
T < T c , whereas for T > T c we expect for the skewness 
(e — > in the thermodynamic limit. It would be inter- 
esting to understand the origins of this behavior of the 
mean-field model analytically. 

For nonzero values of a we find finite-temperate re- 



sults in agreement with the data presented in Sec. IIV CI 
The distributions become Normal in the thermodynamic 
limit. 



V. SUMMARY AND CONCLUSIONS 

We have studied in detail the probability distribu- 
tion function of the ground-state energy of the one- 
dimensional Ising spin chain with random power-law in- 
teractions for several values of the power-law exponent a. 
Using sophisticated parallel-tempering methods (fast for 
small values of a) and a branch, cut, and price algorithm 
(fast for large values of a), relatively large system sizes 
have been studied over the full range of the parameter a. 

For the SK limit, when a = 0, our results agree with 
previous numerical work by Palassini.^ We find by study- 
ing different moments of the distribution, that the SK 
model has a skewed probability distribution function in 
the thermodynamic limit that is well fitted by a modified 
Gumbel distribution, possibly with small Gaussian cor- 
rections. This behavior is not only valid for the ground- 
state energy, but also for energies below the critical tem- 
perature. 

By varying the power-law exponent a we scan sev- 
eral universality classes and show that for the nonmean- 
field regime when a > 0.5 the probability distribution 
functions converge to a Normal distribution in the ther- 
modynamic limit, in agreement with a short-range spin- 
glass ladder. Thus a skewed ground-state energy proba- 
bility distribution function is a characteristic property 
of the mean-field spin-glass model and the change of 
the distribution's characteristic coincides with the tran- 
sition line between the MF and LR universality classes. 
This behavior again poses the question, of whether the 
mean-field description of low-temperature properties of 
spin glasses is adequate for nonmean-field models, as 
has been observed previously by studying other measur- 
able quantities^^^^^i although other studies^ have 
found different results. 

Thus far it is unclear to us why the limiting distribu- 
tion for the SK case is well described by a modified Gum- 
bel distribution with parameter m > 1, i.e., an extreme- 
value distribution for selecting the mth smallest value out 
of a large number of M uncorrelated values^ If all 2 N 
energy levels of a system with N spin were uncorrelated, 
then the ground state would be simply the minimum of 
all 2 N uncorrelated values and a standard Gumbel dis- 
tribution [m = 1) would be the limiting distribution. 
Clearly the energy values of a spin glass are not fully un- 
correlated, but recently it has been observed^ that the 
energy levels of the Edwards- Anderson model behave at 
least locally (i.e., in small intervals) like a random-energy 
model. This might be the underlying reason why a Gum- 
bel distribution seems to describe the data best, as well 
as for the occurrence of a nonvanishing skewness in the 
MF case for a < 0.5. 

In general, we see that by studying the distribu- 
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tions of measurable quantities such as for the ground- 
state energy, we have another approach to discrimi- 
nate mean-field-type behavior from simpler structures 
of the phase space. Therefore this approach supple- 
ments other numerical means of studying the organiza- 
tion of phase space, such as calculating the distributions 
of overlaps^ clustering configurations, 55 or the calcula- 
tion of correlation-matrix eigenvalues; 5 — Hence, it should 
be fruitful to study the distributions of ground-state en- 
ergies in detail also for other models. This is especially 
interesting when a disorder-driven phase transition oc- 
curs, such as for parametrized random bond models, 
random-field systems, or optimization algorithms on ran- 
dom graphs. So far the body of the ground-state energy 
distributions has been tested in detail. More informa- 
tion about the tails of the distributions could be accessed 
using rare-event techniques^ also for the standard spin- 
glass models in finite and infinite dimensions. 
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APPENDIX A: PARALLEL TEMPERING 
GROUND-STATE SEARCH 



In this section we describe the different numerical 
tools to compute ground-state instances of the one- 
dimensional Ising chain fast. As introduced in Refs. 
and |3|] we use parallel tempering Monte Carlo2L2& to 
calculate ground states. In Ref. [l6| it has already been 
mentioned that parallel tempering Monte Carlo performs 
poorly for large values of a. In particular, for a > 2.5 we 
find that in practice the time to find a ground state scales 
exponentially 1 ^ in the system size. In order to overcome 
this limitation we use the branch, cut, and price algo- 
rithm described below. 

In the parallel tempering Monte Carlo method one 
simulates several identical replicas of the system at dif- 
ferent temperatures, and, in addition to the usual local 
moves, one performs global moves in which the tempera- 
tures of two replicas with adjacent temperatures are ex- 
changed. In this way, the temperature of a given replica 
wanders up and down in a random manner, thus provid- 
ing a more efficient sampling of the energy landscape. For 
further details regarding the parallel tempering approach 
see Refs. |3?] and |38[ The parameters of the simulation are 
shown in TableQ] If we take the lowest temperature T mm 
to be 0.05 (T m i n -C T c ), then the minimum-energy state 
found at this temperature is with very high probability 
the ground state. To test whether the true ground state 
has been reached, four criteria have to be met: (i) the 
same minimum-energy state has to be reached from two 
independent replicas at T m i n for all samples, and (ii) this 
state has to be reached during t eq sweeps in both copies, 
(iii) We simulate for further t cq sweeps to ensure that the 
energies found do not change, and (iv) the system has to 
obey the equilibration test for the one-dimensional Ising 
chain, introduced in Ref. In this test the link over- 
lap q\ has to equate the link overlap calculated from the 



FIG. 17: (Color online) Equilibration plot for the one- 
dimensional Ising chain: Average link overlap as a function 
of Monte Carlo steps t eq calculated directly [Eq. HA2I 1. and 
via the internal energy [Eq. IA1I 1 averaged over the last half 
of the sweeps for L — 96, T = 0.05, and a = 0.75. The data 
are equilibrated for t cq « 10 4 MCS, in the simulations 4 x 10 4 
MCS have been used. Data for 2500 disorder realizations. 



internal energy q\(U) via the relation 



91 = 1- 



2T\[U] &V /L\ 



(Al) 



where T C MF is given by Eq. l(T7|). U is given by Eq. 
and 



11- X] (jmfyi [( S i S i) 2 ]^ 
i,j c 



(A2) 



Once both sides of Eq. (|A1() agree, the system is in equi- 
librium (see Fig. I17|) . Note that this is the case for the 
parameters listed in Table [I] If any of the aforemen- 
tioned criteria are not met (usually one instance in 10 ), 
the calculated ground-state instance is rejected. 



APPENDIX B: BRANCH, CUT, AND PRICE 
ALGORITHM 



In this section we briefly explain how exact ground 
states of one-dimensional Ising spin-glass instances can 
be computed fast for large values of a. To this end, we 
extend the branch- and-cut approach to a branch, cut, 
and price (BCP) method originating in combinatorial op- 
timization. Since this approach has not yet been applied 
for spin glasses, we give more details in the following sec- 
tion and discuss the performance in a subsequent section. 
Again, for the fundamentals of the applied algorithm, i.e., 
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TABLE I: Parameters of the parallel tempering Monte Carlo 
simulations. The table shows the total number of Monte Carlo 
steps t cq used for each value of a and L. We use between 
10 and 17 temperatures, depending on the system size, to 
ensure that the acceptance ratios of the parallel tempering 
moves are larger than ~ 0.30. The lowest temperature used 
is 0.05, the highest 1.70. For the internal energy distribu- 
tions (Sec. II V Ell we compute thermally averaged values of 
the internal energy for a given disorder realization after equi- 
librating for t c<1 Monte Carlo steps. The averages are done 
over another period of t eq Monte Carlo steps. For a = 2.50 
and L = 96 the calculations have been done using the BCP 
algorithm (Appendix 0. 
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the standard branch & cut approach, we refer the reader 
to Refs.EOandlil 



1. Algorithm 

The problem of determining a ground state of an Ising 
spin-glass instance is equivalent to determining a max- 
imum cut in the interaction graph associated with the 
system. 40 In the maximum cut problem we are given a 
graph G = (V, E) with nodes V and edges E. The nodes 
correspond to the spin sites, the edges to the bonds. 
Weights dj £ K are given for all edges ij £ E. Let 
W C V be a subset of nodes. The cut 8{W) is defined 
as the set of edges having exactly one endpoint in W. 
The weight of a cut 6(W) is the sum of the weights of 
the edges in the cut, and the maximum cut problem is 
to find a cut 8(W) of G with maximum weight among 
all possible node sets W. Determining a ground state 
of a spin-glass instance amounts to calculating a max- 
imum cut in the interaction graph of the system, with 
edge weights chosen as cy = —Jij for ij £ E. 

The maximum cut problem is NP hard which makes it 
unlikely that there exists a solution algorithm running in 
a number of steps bounded by a polynomial in the size 
of the input. In practice, maximum cuts of reasonably 
sized instances can be determined exactly by using the 
branch-and-cut method from combinatorial optimization 
that has exponential worst-case running time. For an in- 
stance, we always maintain an upper and a lower bound 
for the optimum solution value of the maximum cut. Iter- 
atively we improve upper and lower bounds until they are 
tight enough for proving optimality of a known solution. 



In the upper bound computations, a sequence of linear 
programs is solved. Solving a linear problem amounts to 
optimizing a linear objective function subject to a set of 
linear constraints. Details are explained in Refs. an d 

M 

For an instance of the one-dimensional Ising chain with 
L = 100 spins and a = 3.0, the default version of the 
branch-and-cut algorithm needs roughly 3 h CPU time 
on average on a 1400 MHz Athlon processor. By extend- 
ing the branch-and-cut algorithm to a branch, cut, and 
price algorithm we achieve a better performance. Details 
about pricing algorithms can be found in Ref. 42; 

The underlying idea of a pricing algorithm is as fol- 
lows. There exists a variable for each edge ij £ E, and 
we use the terms edge and variable interchangeably. In 
the pure branch-and-cut algorithm we always work on the 
complete set of variables. However, in the extended algo- 
rithm we start doing branch-and-cut, but only work on 
a small fraction of all variables. We add necessary vari- 
ables (and delete unnecessary ones) dynamically during 
the optimization process. This is done in the so-called 
pricing routine. 

For the one-dimensional Ising chain, we make the as- 
sumption that for big enough values of the parameter a 
the "long-range couplings" between two spins "far apart" 
from each other in the chain do not strongly affect the 
ground state and can be neglected temporarily. Thus, 
we start working on a graph G = (V, E) consisting of all 
nodes but only of a fraction of all edges. In our tests it 
performed best when the input graph consisted of the k% 
edges with highest weights, measured in absolute value, 
where the parameter k is suitably chosen in order to min- 
imize the total running time. (For example, for a = 3.0 
k = 20 is a good choice, for smaller a the value of k is in- 
creased.) At well-defined steps in the algorithm the pric- 
ing routine checks whether there exists a (yet neglected) 
variable that has to be included in the variable set for 
maintaining correctness. If no variable is added, and up- 
per and lower bounds are tight enough, we can prove 
optimality, and stop. For our model, we can further im- 
prove the quality of the upper bound within the BCP 
algorithm by separating not only the cycle inequalities 40 
but also separating heuristically the so-called parachute 
inequalities*^ resulting in an improved bound and an ad- 
ditional speedup. 

When the BCP algorithm is used, solving systems for 
a = 3.0 takes on average 426 ± 55 seconds for L = 100 on 
the same 1400 MHz Athlon processor that needed three 
hours on average for solving the same systems by branch- 
and-cut. In Ref. 1 161 it is reported that parallel tempering 
is less efficient in finding the ground state for bigger val- 
ues of a, because parallel tempering needs longer to relax 
an inconvenient configuration. With the exact algorithm, 
in contrast, we expect pricing to be only effective for big- 
ger a. In this case we expect a speedup by using sparse 
graph techniques as explained above. For small a in- 
stead, the system is of the long-range type, and in the 
worst-case all neglected edges would have to be added in 
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the pricing routine. 

In the following section we experimentally determine 
the running time dependence of the BCP algorithm and 
its dependence on the parameter a and on the system 
size. 



2. Performance of the Algorithm 

In this section we study the performance of the BCP 
algorithm for the one-dimensional Ising chain model. We 
compute ground states of samples for different system 
sizes L and values of the parameter a. The studied ranges 
are a £ {1.0, 1.5, 2.0, 2.5, 3.0}, and L < 96. We compute 
between 1000 and 6000 samples per size and a value for 
small- and medium-sized instances and at least 100 sam- 
ples for the largest instances. All runs are performed 
on a Linux cluster of identical AMD Athlon 1800+ ma- 
chines. Instances of size L < 48 and a > 2.0 are solved 
within seconds; for a — 3.0, computing a ground state of 
L = 280 spins takes on average 5161 ±275 s. The hardest 
instances, L — 96, a — 1.5 needs up to a day computing 
time on one processor. 

As argued before in Ref. IHsl there is no easy and 
"ideal" performance measure for a branch-and-cut algo- 
rithm. This remains true for its extension to the BCP 
algorithm. As a measure of the performance of the latter, 
we could use the needed CPU time which however is ma- 
chine dependent, or the number of solved linear programs 
(lps), see Refs. [5^ and l58l For a < 2 we find that the 
number of lps ni ps is strongly and almost linearly corre- 
lated with the CPU time icpu, see Fig. El The same is 
true for the pure branch-and-cut algorithm. However, for 
a > 2.0 the CPU time for solving a lp considerably varies 
between different samples of the same size, as can be seen 
for L = 64 and a = 3.0 in the scatter plot, Fig. 1181 

A reason for this behavior is the following: In order to 
keep the program flexible, in each iteration we both add 
new constraints to the current linear program and remove 
constraints that once have been added but have turned 
out to be unimportant. (Re-)optimizing a lp is very fast 
if only a small number of constraints changes from one 
iteration to the next but takes considerably longer if a 
substantial change occurs. In the pricing extension, we 
start working on a subset of the variables and might add 
further variables as explained above. Possibly a "bad" 
subset of variables is chosen, in the sense that many of 
the added constraints become unimportant later and are 
removed again. Then the lps change considerably and 
their solution takes long. This is more probable for big 
a, as we start working on a small subset of the variables. 
For smaller values of a instead, we start working on a big- 
ger fraction of all variables and find a stronger correlation 
between number of lps and CPU time«S!l Given the broad 
variation in the CPU time per lp for some values of a, we 
use the mean of the CPU time as a performance measure. 
We notice that the figures remain qualitatively compara- 
ble when the mean of the linear programs is taken instead 
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FIG. 18: (Color online) Scatter plot for the CPU time t C p\j 
in seconds versus the number of lps ni ps for 1000 randomly 
chosen samples with L — 64. The data for a = 2.0 (black 
dots) are strongly correlated. The dashed line is a guide to 
the eye. In contrast, data for a = 3.0 (red crosses) show 
strong sample-to-sample variations. 



of the CPU time. We have checked that the mean of 
both the CPU time and the number of linear programs 
is defined for our sampling as the distribution shows a 
pronounced tail. Performing a detailed statistical analy- 
sis we show that the data are thin-tail distributed^^ with 
a well-defined mean. 

In Fig. E| we show the average CPU time for solving 
an instance as a function of a, for different system sizes 
L. Ground states are computed fast for big values of 
the parameter cr, whereas it takes considerably longer for 
smaller a < 1.5. This effect becomes more apparent with 
increasing system size L. 

We also study the CPU time as a function of the to- 
tal number of edges, i.e., the total number of variables, 
for different values of a. The increase in the CPU time 
with the number of variables is consistent with a poly- 
nomial dependency, even for the smallest studied value 
of cr. When fitting a function of the form f(m) ~ am b , 
with m being the number of variables (bonds), we obtain 
a = 0.009 ± 0.007, b = 1.3 ± 0.1 for a = 2.0. A similar 
behavior can be found when studying the CPU time as 
a function of system size L, see Fig. 1201 

A qualitatively similar behavior can be found in the 
data when plotted as a function of lps instead of CPU 
time (not shown). 
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